<p>
  交易所交易基金（ETF）有许多不同的股票板块和资产类别，为我们提供了各种不同的配对交易选择。我们的数据集包括在纳斯达克或纽约证券交易所进行基金买卖交易的每日数据。
</p>

<p>
  我们使用前3年的数据来选择联系和资产配对的最佳组合（“训练形成期”）。接下来，我们使用2011年至2017年的5年时间（“交易期”）来执行策略。在交易期间，我们使用12个月的滚动窗口数据来获取联系参数（“滚动形成期”）。
</p>

<h3>步骤1：选择配对股票</h3>
<p>
  配对选择的一般方法基于基础分析和统计分析。
</p>

<h4><strong>1)收集可能相关的配对</strong></h4>
<p>
  任何随机对都可以相互关联。这些变量之间可能并没有因果关系，而是由于巧合或存在某种第三种不可见因素而产生了虚假关系。因此，对于我们来说重要的是从一些有共同点的证券开始。为了说明这一点，我们选择了一些在纳斯达克或纽交所交易中流动性最强的基金。这些潜在相关配对之间的关系可能是由于指数、行业或资产类别重叠所造成的。例如，QQQ和XLK就是市场领先指数的两支基金。
</p>

<h4><strong>2)采用统计相关性来筛选交易配对</strong></h4>

<p>
  为了确定分析中所包含的股票配对，将要对预先选择基金对之间的相关性进行分析。以下是我们在统计中常用的三种相关指标：
</p>

<table class="table qc-table">
<thead>
<tr>
<th colspan="3">相关测量技术</th>
</tr>
</thead>
<tbody>
<tr>
<td width="33%">Pearson相关系数</td>
<td>\[r = \frac{\sum (x_i- \bar{x})(y_i- \bar{y})}{\sqrt{\sum (x_i- \bar{x})^2)\sum (y_i- \bar{y})^2)} }\]</td>
</tr>
<tr>
<td>Kendall秩相关系数</td>
<td>\[\tau=\frac{n_c-n_d}{\frac{1}{2}n(n-1)}\]</td>
</tr>
<tr>
<td>Spearman秩相关系数</td>
<td>\[\rho=1-\frac{6\sum d_i^2}{n(n^2-1)}\]</td>
</tr>
<tr>
<td colspan="2"> \(n\) = 每个数据集中的值数
\(n_c\) = 一致数
\(n_d\) = 不一致数
\(d_i\) = \(x_i\)和\(y_i\)对应值秩之间的差异。</td>
</tr>
</tbody>
</table>


<p>
  我们可以在Python中使用SciPy统计值程序库中的函数来获得这些系数。利用训练形成期间的每日股价收益率来计算相关系数。我们可以发现这三种相关性技术对配对基金给出了相同的相关系数排名。Pearson相关性假设两个变量都是正态分布的。因此，本文采用Kendall秩作为相关性测量法，并选择具有最高Kendall秩相关性的配对来进行配对交易。通过使用历史函数并将价格转换为对数收益序列，我们可以得出基金对的每日历史收盘价。设Px和Py为股票x和股票y的历史股价序列，则基金配对的对数收益为：
</p>

\[R_x = ln(\frac{P_{x,t}}{P_{x,t-1}}),   R_y = ln(\frac{P_{y,t}}{P_{y,t-1}})\]   t = 1,2,...,n，其中n是价格数据的数量 

<div class="section-example-container">

<pre class="python">def _pair_selection(self):

    tick_syl =  [["QQQ","XME","TNA","FAS","XLF","EWC","QLD"],
                 ["XLK","EWG","TLT","FAZ","XLU","EWA","QID"]]
    logreturn={}
    for i in range(2):
        syl = [self.AddSecurity(SecurityType.Equity, x, Resolution.Daily).Symbol.Value for x in tick_syl[i]]
        history = self.History(syl, self.lookbackdays,Resolution.Daily)
        # 生成配对股票的对数收益序列
        close = history['close'].unstack(level=0)
        df_logreturn = (np.log(close) - np.log(close.shift(1))).dropna()
        for j in tick_syl[i]:
            logreturn[j] = df_logreturn[j]
    # 不同相关性测量法的估计系数
    tau_coef,pr_coef,sr_coef= [],[],[]
    for i in range(len(tick_syl[i])):
        tik_x, tik_y= logreturn[tick_syl[0][i]], logreturn[tick_syl[1][i]]
        tau_coef.append(kendalltau(tik_x, tik_y)[0])
        pr_coef.append(pearsonr(tik_x, tik_y)[0])
        sr_coef.append(spearmanr(tik_x, tik_y)[0])
    index_max = tau_coef.index(max(tau_coef))
    self.ticker = [tick_syl[0][index_max],tick_syl[1][index_max]]
</pre>
</div>

<h3>步骤2：估算对数收益的边缘分布</h3>
<p>
  为了构造联系，我们需要将对数收益序列Rx和Ry转化为两个均匀分布的值u和v。这可以通过估计Rx和Ry的边际分布函数，并将回归值代入一个分布函数来实现。由于我们对两个对数收益序列的分布没有做任何假设，所以在这里我们会使用经验分布函数来接近边际分布F1(Rx)F1(Rx)和F2(Ry)F2(Ry)。统计模型库中的Python 经验累积分布（ECDF）函数将经验性的累积分布函数（CDF）作为阶梯函数提供给我们。
</p>
<h3>步骤3：估算联系参数</h3>
<p>
  正如上文所讨论的，对于每个阿基米德联系来说，我们通过联系与依附测度Kendall tau之间的关系来估算联系参数。
</p>

<table class="table qc-table">
<thead>
<tr>
<th>联系</th>
<th style="text-align: center;">Kendall's tau</th>
<th style="text-align: center;">参数θ</th>
</tr>
</thead>
<tbody>
<tr>
<td>Clayton联系</td>
<td>\[\frac{\theta}{\theta +2}\]</td>
<td>\[\theta=2\tau(1-\tau)^{-1}\]</td>
</tr>
<tr>
<td>Gumbel联系</td>
<td>\[1-\theta^{-1}\]</td>
<td>\[\theta=(1-\tau)^{-1}\]</td>
</tr>
<tr>
<td>Frank联系</td>
<td>\[1+4[D_1(\theta)-1]/\theta\]</td>
<td>\[arg min\left(\frac{\tau-1}{4}-\frac{D_1(\theta)-1}{\theta}\right)^2\]</td>
</tr>
<tr>
<td colspan="3">\[D_1(\theta)=\frac{1}{\theta}\int_{0}^{\theta}\frac{t}{exp(t)-1}dt \]</td>
</tr>
</tbody>
</table>

<div class="section-example-container">

<pre class="python">def _parameter(self, family, tau):
    if  family == 'clayton':
        return 2*tau/(1-tau)
    elif family == 'frank':
        # debye = quad(被积函数，sys.float_info.epsilon, theta)[0]/theta是一阶Debye函数
    	# frank_fun是平方差
    	# 将frank_fun最小化可以得出直接联系的参数theta
      integrand = lambda t: t/(np.exp(t)-1)
    	frank_fun = lambda theta: ((tau - 1)/4.0  - (quad(integrand, sys.float_info.epsilon, theta)[0]/theta - 1)/theta)**2
      return minimize(frank_fun, 4, method='BFGS', tol=1e-5).x
    elif family == 'gumbel':
        return 1/(1-tau)</pre>
</div>

<h3>步骤4：选择最适合的联系</h3>
<p>
  得到联系函数的参数估计值后，我们使用AIC准则来选择最适合算法初始化的联系。
</p>

\[AIC=-2L(\theta)+2k\]

<p>
  其中\(L(\theta)=\sum_{t=1}^T\log c(u_t,v_t;\theta)\)是对数似然函数，k是参数的数量，在这里k=1。
</p>

<p>
  各联系函数的密度函数如下：
</p>

<table class="table qc-table">
<thead>
<tr>
<th style="text-align: center;">联系</th>
<th style="text-align: center;">密度函数c(u,v;θ)</th>
</tr>
</thead>
<tbody>
<tr>
<td style="text-align: left;">Clayton联系</td>
<td style="text-align: left;">\[(\theta+1)(u^{-\theta}+v^{-\theta}-1)^{-2-1/\theta}u^{-\theta-1}v^{-\theta-1}\]</td>
</tr>
<tr>
<td>Gumbel联系</td>
<td>\[C(u,v;\theta)(uv)^{-1}A^{-2+2/\theta}[(\ln u)(\ln v)]^{\theta -1}[1+(\theta-1)A^{-1/\theta}]\]</td>
</tr>
<tr>
<td style="text-align: left;">Frank联系</td>
<td style="text-align: left;">\[\frac{-\theta(exp(-\theta)-1)(exp(-\theta(u+v)))}{((exp(-\theta u)-1)(exp(-\theta v)-1)+(exp(-\theta)-1))^2}\]</td>
</tr>
<tr>
<td colspan="2">\[A=(-\ln u)^{\theta}+(-\ln v)^{\theta}\]</td>
</tr>
</tbody>
</table>

<div class="section-example-container">

<pre class="python">def _lpdf_copula(self, family, theta, u, v):
    ''' estimate the log probability density function of three kinds of Archimedean copulas '''
    if  family == 'clayton':
        pdf = (theta+1) * ((u**(-theta)+v**(-theta)-1)**(-2-1/theta)) * (u**(-theta-1)*v**(-theta-1))
    elif family == 'frank':
        num = -theta * (np.exp(-theta)-1) * (np.exp(-theta*(u+v)))
        denom = ((np.exp(-theta*u)-1) * (np.exp(-theta*v)-1) + (np.exp(-theta)-1))**2
        pdf = num/denom
    elif family == 'gumbel':
        A = (-np.log(u))**theta + (-np.log(v))**theta
        c = np.exp(-A**(1/theta))
        pdf = c * (u*v)**(-1) * (A**(-2+2/theta)) * ((np.log(u)*np.log(v))**(theta-1)) * (1+(theta-1)*A**(-1/theta))
    return np.log(pdf)
</pre>
</div>

<p>
  所提供最适合的联系应对应于AIC准则的最低值。所选择的配对是“QQQ”和“XLK”。
</p>

<div class="section-example-container">

<pre class="python">self.family = ['clayton', 'frank', 'gumbel']
tau = kendalltau(x, y)[0]  # 估算Kendall秩相关
AIC ={}  # 生成一个关键字为联系系列的字典，值 = [theta, AIC]
for i in self.family:
    lpdf = [self._lpdf_copula(i, self._parameter(i,tau), x, y) for (x, y) in zip(u, v)]
    # 在lpdf列表中用零替换非数字，用有限数目替换无穷大
    lpdf = np.nan_to_num(lpdf)
    loglikelihood = sum(lpdf)
    AIC[i] = [self._parameter(i,tau), -2*loglikelihood + 2]
    # 选择AIC最小的联系
    self.copula = min(AIC.items(), key = lambda x: x[1][1])[0]
</pre>
</div>


<h3>步骤5：生成交易信号</h3>
<p>
  联系函数中包含了两个收益序列依附结构的所有信息。根据Stander Y、Marais D和Botha I，在有联系的交易策略中，使用适合联系推导出\(C(v\mid u)\)和\(C(u\mid v)\)条件边际分布函数的置信区间，即错误定价指数。当市场观察值处于置信区间之外时，表明配对交易机会是可用的。在这里，我们选取95%作为置信区间上边际，5%作为置信区间下边际。置信水平的选择是基于本文的反测试分析，该分析表明，使用95%的置信水平应该可以确定适当的交易机会。
</p>

<p>
  将股票X、股票Y的当前收益设为\(R_x, R_y\)，我们可以将“错误定价指数”定义为：
</p>

\[MI_{X|Y}=P(U\leq u\mid V\leq v)=\frac{\partial C(u,v)}{\partial v}\]

\[MI_{Y|X}=P(V\leq v\mid U\leq u)=\frac{\partial C(u,v)}{\partial u}\]

<p>
  关于进一步的数学证明，请参考Xie W和Wu Y的《以联系为基础的配对交易策略》。通过表1中联系函数的偏导数，可以推导出二元联系的条件概率公式。结果如下：
</p>
<p>
Gumbel联系
</p>
\[C(v\mid u)=C(u,v;\theta)[(-\ln u)^\theta+(-\ln v)^\theta]^{\frac{1-\theta}{\theta}}(-\ln u)^{\theta-1}\frac{1}{u}\]

\[C(u\mid v)=C(u,v;\theta)[(-\ln u)^\theta+(-\ln v)^\theta]^{\frac{1-\theta}{\theta}}(-\ln v)^{\theta-1}\frac{1}{v}\]

<p>
  Clayton联系
</p>

\[C(v\mid u)=u^{-\theta-1}(u^{-\theta}+v^{-\theta}-1)^{-\frac{1}{\theta}-1}\]

\[C(u\mid v)=v^{-\theta-1}(u^{-\theta}+v^{-\theta}-1)^{-\frac{1}{\theta}-1}\]

<p>
  Frank联系
</p>

\[C(v\mid u)=\frac{(exp(-\theta u)-1)(exp(-\theta v)-1)+(exp(-\theta v)-1)}{(exp(-\theta u)-1)(exp(-\theta v)-1)+(exp(-\theta)-1)}  \]

\[C(u\mid v)=\frac{(exp(-\theta u)-1)(exp(-\theta v)-1)+(exp(-\theta u)-1)}{(exp(-\theta u)-1)(exp(-\theta v)-1)+(exp(-\theta)-1)} \]

<p>
  在选择交易配对和最合适的联系后，我们采取下列步骤进行交易。请注意，我们在每个月的第一天使用过去12个月的每日数据来实施步骤1、2、3、4，这意味着我们的经验分布函数和联系参数theta估值会每月更新一次。在每月总结中：
</p>

<ul>
 	<li>在12个月的滚动形成期内，使用日收盘价计算这对基金的每日对数收益，然后计算Kendall秩的相关性。</li>
 	<li>估算X和Y对数收益的边际分布函数，分别为ecdf_x和ecdf_y。</li>
 	<li>将Kendall tau代入联系参数估计函数，得到theta的值。</li>
 	<li>对两个价格系列进行线性回归。这个系数用来决定股票X和股票Y的买卖数量。例如，如果系数是2，每买入或卖出一股X，就买入或卖出2股Y。</li>
</ul>

<div class="section-example-container">

<pre class="python">def _set_signal(self):
    history = self.History(self.ticker, self.lookbackdays,Resolution.Daily)
    # 生成配对股票的对数收益序列
    close = history['close'].unstack(level=0)
    logreturn = (np.log(close) - np.log(close.shift(1))).dropna()
    x, y = logreturn[self.ticker[0]], logreturn[self.ticker[1]]
    # 估算每个交易日的Kendall秩相关
    tau = kendalltau(x, y)[0]
    # 估算联系参数：theta
    self.theta = self._parameter(self.copula, tau)
    # 模拟两对股票收益的经验分布函数
    self.ecdf_x, self.ecdf_y  = ECDF(x), ECDF(y)
    # 对两个历史收益序列进行线性回归
    self.coef = stats.linregress(x,y).slope
</pre>
</div>

<p>
  最后，在交易期间，我们每天使用经验分布函数ecdf_x和ecdf_y将当天的收益转换为u和v。然后，利用估算的联系C，在每个交易日计算两个错误定价指标。\(MI_{Y|X}&lt;0.05\)和\(MI_{X|Y}&gt;0.95\)的日子，算法构建了X的空头头寸和Y的多头头寸。\(MI_{Y|X}&gt;0.95\)且\(MI_{X|Y}&lt;0.05\)时，构建Y的空头头寸和X的多头头寸。
</p>

<div class="section-example-container">

<pre class="python">def OnData(self,data):
    for i in self.syl:
        self.price_list[i].append(self.Portfolio[i].Price)
    # 计算两支股票的当天对数收益
    if len(self.price_list[self.syl[0]]) &lt; 2 or len(self.price_list[self.syl[1]]) &lt; 2: return
    else:
        return_x = np.log(float(self.price_list[self.syl[0]][-1]/self.price_list[self.syl[0]][-2]))
        return_y = np.log(float(self.price_list[self.syl[1]][-1]/self.price_list[self.syl[1]][-2]))
    # 使用经验分布函数将两种收益转换为统一值u和v
    u_value = self.ecdf_x(return_x)
    v_value = self.ecdf_y(return_y)
    # 用估算联系来计算u和v的错误定价指数
    self._misprice_index(self.copula, self.theta, u_value, v_value)
    quantity = self.CalculateOrderQuantity(self.syl[1],0.4)
    if self.MI_u_v &lt; self.floor_CL and self.MI_v_u &gt; self.cap_CL:
        if self.Portfolio[self.syl[0]].Quantity &lt; 0 and self.Portfolio[self.syl[1]].Quantity &gt; 0:
            self.Liquidate(self.syl[0])
            self.Liquidate(self.syl[1])
            quantity = self.CalculateOrderQuantity(self.syl[1],0.4)
            self.Sell(self.syl[1], 1 * quantity )
            self.Buy(self.syl[0], self.coef * quantity)
        else:
            self.Sell(self.syl[1], 1 * quantity )
            self.Buy(self.syl[0], self.coef * quantity)
    elif self.MI_u_v &gt; self.cap_CL and self.MI_v_u &lt; self.floor_CL:
        if self.Portfolio[self.syl[0]].Quantity &gt; 0 and self.Portfolio[self.syl[1]].Quantity &lt; 0:
            self.Liquidate(self.syl[0])
            self.Liquidate(self.syl[1])
            quantity = self.CalculateOrderQuantity(self.syl[1],0.4)
            self.Buy(self.syl[1], 1 * quantity )
            self.Sell(self.syl[0], self.coef * quantity)
        else:
            self.Buy(self.syl[1], 1 * quantity )
            self.Sell(self.syl[0], self.coef * quantity)
</pre>
</div>
